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Abstract 

We introduce a nonstationary spatio-temporal model for gridded data on the sphere. The 
model specifies a computationally convenient covariance structure that depends on hetero¬ 
geneous geography. Widely used statistical models on a spherical domain are nonstationary 
for different latitudes, but stationary at the same latitude (axial symmetry). This assump¬ 
tion has been acknowledged to be too restrictive for quantities such as surface temperature, 
whose statistical behavior is influenced by large scale geographical descriptors such as land 
and ocean. We propose an evolutionary spectrum approach that is able to account for 
different regimes across the Earth’s geography, and results in a more general and flexible 
class of models that vastly outperforms axially symmetric models and captures longitudinal 
patterns that would otherwise be assumed constant. The model can be estimated with a 
multi-step conditional likelihood approximation that preserves the nonstationary features 
while allowing for easily distributed computations: we show how the model can be fit to 
more than 20 million data points in less than one day on a state-of-the-art workstation. The 
resulting estimates from the statistical model can be regarded as a synthetic description (i.e. 
a compression) of the space-time characteristics of an entire initial condition ensemble. 
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1 Introduction 


Providing efficient and flexible models for data on a spherical domain is a topic of great 
importance in climate science, as the statistical model can be used to fit global data. In par¬ 
ticular, in the context of Earth System Models (ESMs), this could lead to efficient methods 
for compressing large quantities of data. Isotropic mo dels have bee n widely acknowledged 


as being inadequate for data on a spherical domain (IGneitind . 


2013af ). and defining valid 


nonstat ionary pr ocesses is listed among the sixteen open problems in modeling spherical 


data in 


ferential equation (ILindgren et al. 


Gneitind (120 13 m. By regarding a random field as solutio n of a stochastic par 


2 0111 ) on a spherical domain, 


Bolin and Lindgrcn 


ial dif- 


( 20111 ) 


proposed a nested stochastic partial differential equations approach, which yielded a field 
with Matern-like covariance structure but could also be extended beyond axial symmetry 
to no nstationary models by allow i ng flexible differ e ntial o perat ors in the stochastic equa¬ 


tion. 


Jun and Stein ( 2007 


20081 ): 


Jun et al. 


( 120081 ): 


Jun (1201 if ) restrict three-dimensional 


isotropic fields to a sphere and apply partial derivatives with respect to latitude and lon¬ 


gitude, obtaining a model which assumes stationarity if th e dat a are at the same latitude, 

a), 


and n onstationarity 


H i tczenko and Stein ( 2012 ); 


otherw i se (axially symmetr ic dJonesl . 


Huan g et al. 


1963b see theoretical details in 


020121 ) b Such models are conceptually attractive 


f or data su ch as surfa ce te mpera 


Castruccio and Stein (120131) and 


u re, w h ose statis tical pr opert ies clearly depend on latitude. 


Castruccio and Gentonl ( 2014 ) proposed a spectral approach 


that is flexible and computationally efficient when the data are on a regular grid over the 
sphere. This method proposes to separately consider the process by latitudinal bands, fit 
a Matern-like covariance across longitudes, and then estimate the multi-band dependence, 
thus reducing the likelihood evaluation with the full dataset to a low dimensional param¬ 
eter space. The main limitation of these models, as acknowledged in the aforementioned 
literature, is the assumption of stationarity in longitude at each latitude. 
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For physical quantities such as surface temperature, it is expected that large scale ge¬ 
ographical descriptors such as land/ocean will impact the statistical behavior of the data. 


Recently 


■Tunl ( 20141) proposed a modified Matern process with smoothness changing over 


land and ocean, which showed dramatic improvements over the axially symmetric model. 
The model parameters, however, were not simple to interpret given their definition through 
a differential operator over an isotropic process and the fitting procedure was not computa¬ 
tionally feasible for analyzing millions of data points. 

This work introduces a new class of covariance functions on spheres that includes axially 
symmetric models as special cases and is capable of incorporating geographic covariates 
into the model. For the surface temperature data we consider, the most prominent and 
influential large scale geographic descriptor is land versus ocean, so we focus our work on 
this covariate, but as we describe in Section 13.21 the ideas generalize to other covariates as 
well. We also propose a reformulation of the latitudinal dependence of the model in terms 
of a stationary AR(1) process, and introduce a nonstationary generalization which is more 
flexible in capturing different behaviors in the tropics. 

For inference, we devise a step-wise conditional likelihood approach that fully exploits 
the gridded geometry of climate model output and is able to achieve a fit of more than 20 
million data points in less than one day by allowing code parallelization on a state-of-the-art 
workstation. The proposed method vastly outperforms the axially symmetric model in terms 
of standard model selection metrics, and is also able to capture patterns in the longitudinal 
contrasts that would be otherwise assumed constant. The set of estimated parameters can 
then be used to almost instantaneously produce new surrogate simulations on a common 
laptop, thus allowing an end user to conveniently test initial scientific hypotheses on a high 
(spatial) resolution ensemble without downloading it, or remotely aggregating data in space 
or time and losing valuable information at fine scale. Besides, the estimated parameters 
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can be re garded as descriptors o f the i nformation for every member of the gi ven initial 


19891; 


ensemble flCastruccio and Genton . 


Hansen and Yul . 


20161 ). and thus as a compression algorithm (IRissanenl . 


2001J). The proposed statistical model achieves a compression rate of 


approximately 3:100, which is vastly superior to traditional bit-by-bit compression algorithms 
that can achieve at most a 1:5 ratio. 


The model can also be viewed as an emulator of an initial condition ensemble (ICastruccio and Genton . 


2016b under t 


(ILqrenzJ, 


1963 


re assumption that run s are indepen d ent for different init i al con ditions 


Collins and Allen, 


2002 


Collins, 


2002 


B ranstator and Tengl, 


2010) • The 


use of emulators as data compressors is, to our knowledge, new to the cli mate co nnnu- 


scenario extrapolation flHolden and Edwards! . 


2010 


Holden ct al. 


2013 


nitv as thev are tradit 

ionallv used for calibra 

■ion and sensi 

tivitv analvsis fSanso et al.. 

OO 

0 

0 

CM 

Sanso and Forest. 

20 oJ: 

Bhat et ah. 

2012; 

Drignei et ah. 

2008; 

Chang et al.. 

2015 

) or 


Castruccio et al. 


2 0141 ). Having statistical models (emulators) that accurately describe the model output 


allows us to avoid storing the entire initial condition ensemble, whose individual member 
requires significant storage space. This proposed methodology has shown promising results 
and can be generalized to multiple variables (not necessarily in the atmospheric part of 
the model), climate models and scenarios, and to finer temporal scales. In all these cases, 
the benefits of a statistical-based data compression will be even more evident as the size 
of the data, and consequently the expected time of downloading the full climate run, will 
significantly increase. 

The remainder of the paper is organized as follows. Section [2] introduces the data set 
and discusses nonstationarity across longitude. Section [3] describes the statistical model, 
discusses the computational challenges that arise when fitting this with very large data 
sets and suggests a stepwise model-fitting approach to address these challenges in a way 
that exploits the geometry of the sphere. Section [4] shows the comparison with the axially 
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symmetric model. Section [5] shows how the fitted model can be used to compress the initial 
condition ensemble and how to generate surrogate runs from the estimated parameters. 
Section EH concludes with a discussion. 


2 The CMIP5-CCSM4 ensemble 


The Coupled Model Intercomparison Project phase 5 (CMIP5 


Ta ylor et al. 


2012) is a set 


of coordinated experiments from many modeling groups to provide uniform and comparable 
assessment of climate response under different climate models for the latest IPCC Assessment 


Report (IPCC 


Model 4 (CCSM4 


(rcp85 


20 m In this wor 


Gent et al 


Van Vuuren et al 


c we focus on the NCAR Community Climate System 


2 0111) . under the Representative Concentration Pathway 8.5 


2 0111 ) from 2006 to 2100, for a total of 95 years. We consider 


annual temperature at surface (considered at a standard height of 2 meters above ground 
level), which is on a regular 192 x 288 grid over latitude and longitude. We remove the 
bands near the poles (south of 62 degrees south and north of 70 degrees north) so that 
each spatial field consists of 142 x 288 points. The removal of the Arctic and Antarctic 
bands was performed to avoid inferenc e on latitudes where two loc at ions at one longitudinal 


l ag we re ve ry c l ose, con s istent l y wi th 


(20 141 ) and 


Castruccio and Stein (120131) : 


Castruccio and Genton 


Castruccio and Gentonl (j2016|). This would have led to very smooth spatial 


processes, and consequently computational challenges that would have added to the already 
substantial complexity of the inference scheme. Under rcp85, the CCSM4 was run un der 6 


different sets of initial conditions, the r efore generating 6 independent re alizations (ILorenzI . 


1963 


Collins and Allen. 


2002 


Collins, 


2002 


Br anstator and Teng. 


2010 ). The total size of 


the dataset is therefore 142x288x95x6 = 23.3 million points. The movie movie_cm.avi in the 
supplementary material shows a realization of this climate model. The annual temperature 
shows evidence of normality, as reported in the supplementary material. 
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We denote by T r the temperature for realization r = 1,..., R, by L m e (—vr/2, 7 t/2 ),m = 
1 ,M the latitude, by E n = 2irn/N, n — 1,..., N the longitude, by tk,k = 1,..., K the 
year, where R — 6, M = 142, IV = 288 and K = 95. Thus, the temperature for realization r 
is represented as 


T r — {T r (Li, Ei, ti ),..., T,.(Lm, E\, ti), T r (Li, t 2 , ti),..., T t {Lm, En , £&:)}• 


The defining assumption of axially symmetric models is that the process is stationary 
across longitude at each latitude. In this work, we relax this assumption to allow geographic 
covariates to be incorporated into the covariance function and inform more complex spatial 
dependence structures. Local geography can have a strong impact on the statistical charac¬ 
teristics of surface temperature data, so a natural deviation from the stationary assumption 
is to allow the statistical properties of the process to differ over land and ocean, which is the 
most dramatic geographic descriptor at large scales. A simple modeling solution is to divide 
the temperature at each location (L m ,E n ) by the standard deviation SL m ,f n obtained from a 
simulation from the same climate model but wit h no forcing (a control run in geoscience ter¬ 


minology), as proposed in 


Castruccio and Steiru (120131 1. to obtain more realistic conditional 


simulations. While producing improved results, this approach does not allow for a chang¬ 
ing correlation structure across longitude, and in particular across land and ocean. Indeed, 
empirical estimates of the second-order (covariance) structure show a strong dependence on 
the land/ocean variable. To see this, we consider the difference between two realizations 
Ti — T 2 (to remove any trend), normalize it by SL m ,i n , and compute 


k =1 


N 


5> j (4 


Ti(L m , E n , tk) — T 2 (L m , E n , tk) _j 


( 1 ) 


n= 1 S Lm/n 

for j = 1,2, which is the periodogram of a tapered version of the data averaged over time at 
latitude L m , where the taper h 1 is a smooth function that is equal to zero when E n 6 ocean, 
and h 2 is a smooth function that is equal to zero when E n G land, and p(h J ) is a normalizing 
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constant. Thus we can view 1 2 as a periodogram for the land data averaged over time at 
latitude L m , and |/| 2 as a periodogram for the ocean data averaged over time at latitude 
L m . In Figured] we plot log(| fl '/') at latitude L m = 41°. Because the two log periodograms 
in Figured! are not parallel-which would indicate similar correlation structure over land and 
ocean-it is clear that the data exhibit land/ocean nonstationary correlation, and that the 
process over the ocean is much smoother than the process over the land at this latitude. 


land/ocean mask across longitude 



Figure 1: Comparison of the land/ocean periodogram of the difference between two realiza¬ 
tions of rcp85 computed with (dj), each pixel being normalized by its standard deviation from 
the control run. The latitude band represented is ~ 41°N, and the periodogram is averaged 
over all years. On the top the land/ocean mask across longitude is represented. 

The land and ocean periodograms in Figure [D motivate the use of a model that allows 
for a different behavior across land and ocean for the global space-time temperature data. 
In this work, we define a nonstationary model with changing spectrum across these two 
domains, denoted as evolutionary spectrum. The details are provided in Section 13.21 
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3 The space time model 


In this section we describe the global space-time model. We first introduce in Section 13.11 
a fundamental result that allows us to fit the stochastic component of the statistical model 
without defining a parametric form for the mean. In Section [3T2l we describe the model. Sec¬ 
tion 13.31 shows the results and discusses the computational challenges of fitting the proposed 
model on a dataset with tens of millions of data points. 

3.1 Preliminaries 


Denote by E(T r ) = /x the mean temperature across realizations. Since realizations differ 


just in their initial condition, and sin ce climate models tend to forg e t the i r initia 


after a short number o: 


Branstator and Tengj . 


temp oral steps flLorenzI . 


1963; 


Collins and Allen. 


2002 


Collins 


state 


2002 


2010l ). we can assume that the space-time field T r is independent 


across r: 


T r = n + e r , e r ~ Af(0, £(0)), 


( 2 ) 


where 6 is a vector of unknown covariance parameters. The noticeable advantage of hav¬ 


ing independent realizations is t hat 6 can be est i matec 


via restricted loglikclihood. 


without any parametrization of pi 


Castruccio and Stein (120 130 proved the following result, which 


formulates the restricted loglikclihood for T r in a computationally convenient form. 


Result 1 Denote with A T r = T the average temperature across realizations. Let 

D, r = T r — T. The negative restricted loglikelihood for ((2D is 

((0;D) = ™ J f- 1 > fo g ( 2n) + i(R - l)fog|S(g)| 

+i KNMlog(R) - i I), 1 >:(#) '!), 

Also, the corresponding estimator for /j, obtained by generalized least squares is ft = T. 


Throughout this work we make use of ([2D to estimate the space/time structure of the data. 
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3.2 Sphere-Time Covariance 


Denote by e(4;r) = {e r (Li, l\, t k ), ..., e r (Ljv, Em, 4)} the vector of the stochastic term of 
(|2|) at time 4, by T (t k ',r) the temperature at year tk for realization r and by D (t k ',r) = 
T (tk'-, r) — T (t k\ r). We assume that £(t k ',r ) is correlated across time, and previous work 


(jCastruccio and Gentonl . 


20141 ) has shown that an AR(2) model with different coefficients 


for every grid point is sufficiently flexible, as no evidence of cross-temporal dependence or 
nonseparability between space and time was found on annual scale: 


e(4; r) = $1 e(4 - 1; r) + $ 2 e(4 - 2; r) + SH(t fc ; r), 


(4) 


where and $ 2 are two NM x NM diagonal matrices with the autoregressive coefficients 
for each location, and S is a NM x NM diagonal matrix with the associated standard 
deviations. 

The unsealed innovations H r (L m ,£ n ,t k ) are independent across time and describe the 
spatial dependence across the sphere. We propose to model the process in the spectral 
domain across longitudes, and then to model the dependence across latitudes: 

>JV-1 


H r(L m J n ,t k ) = J2c=o /Wn( c )e nC Hr(c,Tn,4), 

corr |h r (c,L m ,t k ), H r /(c', L m >, f fe /) j = l{c = c',k = k',r = r'}p Lm ,L m ,(c), 

where c is a wavenumber. If = /z, m f° r every L m , than this would be a standard sta¬ 

tionary model with spectral density |/i m | 2 . This proposed model assumes that the spectral 
density is not exactly constant in longitude, but evolves (hence the term evolutionary spec¬ 
trum) according to |/L m ,£ n | 2 - Models with evolutionary spectra allow us to flexibly specify 
the local covariance properties at every location while ensuring that the resultin g cov a riance 


function is positive definite. Evolutionary sp ectra were first introduced by 
to model nonstationary time series data, and 


Guinne ss and Stein (2013j) provided computa- 


Pr iestl evl ( 19651 ) 


tionally efficient methods for fitting models with evolutionary spectra to nonstationary time 
















series data. In this work, we adapt the evolutionary spectra to model nonstationarity over 
longitude rather than over time, which requires a discrete spectrum because of the circular 
domain on which the data at a single latitude fall. 


Certain regularity conditions can be imposed on ./b m / n near the poles to achieve mean 
square continuity (see supplementary material). PL m ,L m , is the correlation in the spectral 
domain (or coherence in spectral analysis) between latitudes L m and L m / among H r for 
the same wavenumber, time and realization. Alternatively, one could deviate from the 
stationary assumption across longitude by introducing dependence across wavenumber in 
H r . Our choice to use evolutionary spectra to model the nonstationarities is motivated by 
the need to include geographic cova riates. 


The typical approach flPri estlevl . 


19651 ) for fLm,e n {c) is to describe the dependence on 


according to covariates X^(L m ,£ n ) as 


/w.M = £/L(c)^(£m,U 

j = 1 


( 6 ) 


In this work, we propose a novel model where land and ocean are included as covariates 
to allow for different statistical behaviors across these two domains. Thus, fL m ,e n (c ) can be 
expressed as 


h m ,e n (c) = fl m (c)b land (L m ,£ n ) + /l m (c) {l - 5 land (T m ,£„)} , 


(7) 


where the component spectra are modeled according to the parametric form 


l/L(c)l 2 = 


i> J L 

J-'r 


{K.) ! +4 


S mu ] Hr i " +i/2 ’ 


J = 1,2, 


( 8 ) 


and 6| and is a function between 0 and 1 that modulates the relative contribution of the land 
regime. (J8J) is a Matern-likc spectrum, which is modified for the case of data on a circle 
to allow for a smooth transition at high wavenumbers and has been shown to adequately 
capture the longitudinal behavior of temperature at surface for different latitudes better 
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than the traditional Matern-like spectrum flCastruccio and Steinl . 


2013 


2014). 


Poppick and Ste in]. 


Choosing an indicator function for 6 j an( j would result in abrupt transitions between the 
two regimes at land/ocean boundaries and in misfit of the data, as shown in the supplemen¬ 
tary material. We therefore introduce a smoother taper function to transition between land 
and ocean: 


• Let /,„(£„) denote the indicator function of land at latitude L m and longitude i n . 
Wherever there is a land/ocean transition, we modify I m (£ n ) so that is equal to one 
for gi m more grid points, where gi, m is an integer number that can also be negative. 
The modified indicator is denoted by I m (£ n ] g Lrn ). 


Compute the Tukey taper function (Tukev, 


19671) with range 7 L m - 


1 + cos 1 ^-( £ 


7l„ 


7i ra ) < 


= 1 


1 + cos {- 1 


7 L„ 


0 < 7 L m < , 

7 n m / 2 < 4 < 1 - 7Um/ 2 , 
, l-7i m /2<4<27r. 


(9) 


• Convolve Im{^n',gL m ) witl1 Wm(4;7i m ) : 

N 

^land^ m ’ ^ n ’ 9L m , 7 L m ) = /* ] 4(4; gL T ^)w m {£ n — 4'j 7 h m )- (10) 

n'=l 

This formulation imposes a symmetric land/ocean transition (i.e. land/ocean and ocean/land 
transitions are equally smooth); however, more sophisticated models with asymmetric tran¬ 
sitions have been tested but have not yielded significantly better results. Similarly, no 
significant improvements have been observed if a different taper is used (as shown in the 
supplementary material) or and 7 L m are assumed different across oceans. If we con¬ 
strain 



1 ot r — ot t 


Vt = Vj 


( 11 ) 
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then in (ED f Lml „ = ./> and t 


Castruccio and Stein! (120131) : 


le model becomes stat iona ry across longitude. 


Castruccio and Gentonl (20141. 


20161 ) propose the following 


parametric model for pL m ,L m , ( c ) hi (jSJ): 
PLm,L m , (c) = PL m -L m , (c) = 


\L m — L m , I 


= <p(c) 


I Lm Lm' | 


( 12 ) 


{l+4 S in 2 (+,r)} 

with ip(c) = | 1 + 4 sin 2 ^ c n)V ' This process is equivalent to the following AR(1) process in 


latitude: 


H Uc) = 


e Lr, 


<^(c)H L m _!(c)+ e Lm (c), m = 2, 
e Ll (c)~J\f( 0,1), m = 1, 

~ Af(0,l — (p(c) 2 ), m> 1 , 


(13) 


where var(eL m (c)) = 1 — ’p(c) 2 to allow unit variance on Hz, m (c). While the coherence in 
(lH has been previously used in literature, the formulation of the latitudinal dependence in 
terms of an autoregressive process has never been acknowledged. 

The formulation of the dependence as a stationary AR(1) process allows for generalization 
to nonstationary latitudinal processes to increase the model flexibility. In particular, in 
addition (I12p we also propose a novel nonstationary AR(1) model for the coherences, with 
latitudinally-varying autoregressive parameters, that is 


<PL m {c) = 




I T Lr, 


(14) 


{1+4 sin 2 (£*)}’ 

Our diagnostics have shown that the coherences are nearly stationary outside of the tropics, 
so we fit nonstationary coherences within —23° < L < 23° (i.e. in the tropics), while we 
assume a constant outside this region. 

Thus, the model consists of three sets of parameters to be estimated 

• The temporal parameters, consisting of all the entries in 3>i, $2 and S in (ED, which 
will be denoted as 0 t ime- 


The longitudinal parameters, consisting of , a 3 Lm , z/) jm; 


111 


and (c/L m ,7nJ in 


(ROD for m — 1,..., M. We denote the collection of all parameters as 0i on - 
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• The latitudinal parameters, consisting of (£_L m , ry m ) in (fT4|) m = 1 ,..., M. We denote 
them as 0i at . 

3.3 Model fit and computational considerations 

Despite the computationally convenient form in (jHJ) (which can be further simplified as shown 
in the supplementary material) it is not feasible to perform a global optimization, since this 
would imply maximizing the likelihood over more than 100,000 parameters (the temporal 
part requires 3 parameters for each of the 142 x 288 ~ 41, 000 locations, the spatial part 
requires a total number of parameters shown in the first row of Table [T]). We therefore propose 
successive conditional approximations of (]3]) by assuming independence across increasingly 
large subsets, each approximation assuming the parameters from previous steps to be known 
and fixed. 

1. Estimate the temporally autoregressive parameters 0 t ime> assuming that the innova¬ 
tions H(4; r) are independent across latitude and longitude. 

2. Consider 0 time fixed at their estimated values and estimate 0i on , assuming the innova¬ 
tions H(t fc ; r) are independent across latitudes. 

3. Consider 0 t ime and 6\ on fixed at their estimated values and estimate 0i at . 

The choice of the blocks in the approximation, as well as the approximation order is dic¬ 
tated by the geometry of the problem as well as from physical considerations. Estimating the 
temporal structure for each location assuming no cross-correlation allows for a very fast (and 
scalable) computation of approximation 1. The choice of latitudinal bands in approximation 
2 allows flexible estimation of the statistical parameters across latitude, which is the main 
descriptor of surface temperature. Further, this choice results in exactly circulant matrices 
across longitudes in the axially symmetric case, a feature that allows very fast computa¬ 
tions in the spectral domain. This conditional approximation scheme can be generalized to 
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allow for vertical profile of temperatures as in 


Castruccio and Gent on (120161 ). and can also 


be applied to any large space-time data set where the geometry, as well as the physics of 
the problem su ggest t he bl ocks, e.g. regions of interest in functional Magnetic Resonance 


Imaging (ICastruccio et al 


20161) 


The sequential model-fitting procedure can also be used to fit axially symmetric versions 
of the model for the innovations. This involves imposing the constraint (TTT1) in approxima¬ 
tion 2. In Figure [2] we see a comparison of the evolutionary spectrum model with the axially 
symmetric model (i.e. with constraint (ITT]) ) in terms of estimated parameters and loglike- 
liliood. Figure [2ja-c shows how land and ocean parameter estimates for the evolutionary 
spectrum are very different from those of the stationary model, and how there is a consistent 
difference across latitude. Figured shows how the smoothness parameter is smaller for land 
than for ocean which implies, as noticed in Figured] that ocean temperatures tend to have a 
smoother behavior across the same band compared to land. Given the very large size of the 
data set, the parameter estimates are very precise and the estimated standard deviation^ are 
two orders of magnitude smaller than the point estimates. Thus, we chose not to report the 
confidence intervals as they are very small compared to the differences across latitudes. We 
also report in Figure [2tl the individual maximum loglikelihoods for each band. The loglikc- 
lihood shows a noticeable improvement for the evolutionary spectrum approach, especially 
in the Southern Hemisphere. In latitudes where there is no land, such as the southernmost 
bands considered (we removed the Antarctic regions) the evolutionary spectrum and the 
axially symmetric model are the same and thus have the same loglikclihood. 

Approximation 3 would require estimating £ Lm and when —23° < L < 23° and a 
constant value for both parameters outside the tropics, for a total of 50 parameters. Since 
a likelihood maximization for such number of parameters and with tens of millions of data 


Computed at each step conditional to the previous steps. For example, the standard deviations for 
(</>£, a J L , v J L ) are estimated conditional on the temporal parameters. 
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latitude 



Figure 2: Comparison of the models with evolutionary spectrum (J7J) and the axially sym¬ 
metric model with constraint (1111) in terms of (a) logand log(0^ m ) for j = 1,2, (b) 
a Lm and oP L ■, (c) z>£ m and and (d) loglikelihood. A smoothing spline has been applied 
to the estimated parameters for the evolutionary spectrum approach in a-c since the pattern 
were less regular. 


points is not feasible, we first consider (fT4j) for adjacent bands, obtain their estimates in- 
dependency for every pair of bands, which we denote as and , and consider these 
estimates as fixed in approximation 3. (Since every band is involved in two fits, by conven¬ 
tion at latitude L m we plug in the estimates from bands ( L m , L m+1 )). The fitted parameters 

A/n\ /o\ 

for (TT3j) and (fT4j) . along with Q jm and t) J, can be found in Figure S4 in the supplement. 
The stationary model shows some misfit, especially for this is due to model assuming 
a constant value across latitude for the coherence, while this is significantly smaller in the 
southernmost regions and at some tropical latitudes. The parameters of the nonstationary 
AR(1) model (TT4l) instead are fixed and equal to an d t l ' in the tropical regions (by con¬ 
struction) and, while still not capturing nontrivial latitudinal patterns outside the tropics, 
it results in a larger and more satisfactory estimate for £. 
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4 Model Comparison 


Tabic Q] shows a comparison among a model that assumes spatial independence (denoted 
ind ), the axially symmetric model (denoted ax ), a model with land/ocean evolutionary 
spectrum with a stationary AR(1) latitudinal process (TT2l) (denoted ev-st ) and one with a 
nonstationary latitudinal AR(1) process (TT4|) (denoted ev-nst ). 


Table 1: Comparison between different models in terms of number of parameters (excluding the 
temporal ones), compu tational time (hours), normalized restricted loglikelihood ©, and Bayesian 
Information Criterion ( Schwarj . 1978 1). 


Model 

ind 

ax 

ev-st 

ev-nst 

4 pararn 

0 

428 

1138 

1234 

time (hours) 

1.4 

1.5 

13.8 

14.8 

Aloglik/NMT(R-l) 

-2.87 

-0.61 

-0.0018 

0 

BICxlO 8 

-0.1677 

-1.0465 

-1.2832 

-1.2839 


The model assuming independence is clearly the fastest to fit, as once the temporal part 
is estimated, the full likelihood can be evaluated just once. The axially symmetric model 
requires spatial parameters, but the computational time is almost equivalent and the im¬ 
provement both in terms of normalized likelihood and BIC is noticeable. The evolutionary 
spectrum model requires approximately three times more parameters than the axially sym¬ 
metric model and a noticeable increase of computational time (mostly because of the 2-band 
step). The resulting model, however, shows a dramatic improvement; the loglikelihood im¬ 
proves by 0.6 units per observation, and improves the BIC despite the large increase in the 
number of parameters, ev-nst requires more parameters (the plug-in estimates and z/j 2 ^ 
at the equator) and there is small indication of a further improvement in the fit. 

To assess the quality of the fit, we compute the contrast variance 


A, 


m EfcLi Eti { Kr(L m , 4 , 4 ) - h r {L m , 4-i, 4 )} 2 , 


ew;m,n 

^ns;m,n ~KR ^ y f-— i ^ v -i —i | 


(15) 
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Figure 3: Estimated and fitted variances for the contiguous differences at approximately 28° 
North for different longitudes, averaged across time and realizations, (a): A ew;m! . and (b): 
A ns;m) . as in (fT51) . The vertical axis is plotted on a log scale. 

where ew=east-west and ns=north-south, and compare them with the corresponding fitted 
quantities according of ax (axially symmetric) and ev-nst (latitudinally nonstationary evo¬ 
lutionary spectrum). The result for a chosen band at approximately 28° North is shown 
in Figure [3J Both panels show the limits of axial symmetry which, assuming longitudi¬ 
nal stationarity, results in a constant value across longitude. This is clearly not adequate 
for temperature data, as there are significant longitudinal patterns generated by different 
land/ocean domains; for this latitude, A ew;mj . is nearly ten times larger over land than it 
is over ocean. The evolutionary spectrum model proposed here is noticeably more flexible 
and is able to accurately capture the changes across longitude in the contrast variances. It 
is apparent how different domains have different behaviors, and thus different spatial corre¬ 
lation, and how the fitted ev-nst allow for a smoother spatial behavior over the ocean. The 
evolutionary spectrum proposed is particularly effective in capturing A ew; m,- hi Figure [3^, 
while some misfit is present in the Pacific Ocean for the north-south contrast variances in 
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Figure [3b. 


longitude latitude 

-150 -100 -50 0 50 100 150 -60 -40 -20 0 20 40 60 



Figure 4: Estimated and fitted averaged contrasts variances as in (fT6lh (a): A^.. (b): A^.., 
(c): A}™. and (d): A^... The vertical axis is plotted on a log scale. 


(16) 


To assess the fit for multiple latitudes, we compute the average contrast variance 

A lat = A- 

Alon _ J_ a 

M A= 1 L -^j\m,ni 

where j = (ew, ns}. In Figure [Ur-b, the values for j = {ew} contrasts are shown, and while 
both ax and ev-nst are able to capture longitudinally averaged variances in panel (b) (apart 
from a misfit of both models in the southernmost bands), only ev-nst is able to capture the 
pattern in latitudinally averaged variance in (a), since ax assumes constant variance across 
longitudes. Figure 0]>d shows the values for j = {ns}. Similar remarks as in the previous 
case hold, but ev-nst does not fully capture the patterns of latitudinally averaged variance 
in panel (c), while the two models performs similarly (and with some degree of misfit in the 
southernmost latitudes) in longitudinally averaged variances in panel (d). 
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5 Simulating the initial condition ensemble 


We now proceed with simulating surrogate (emulated) runs according to the evolution¬ 
ary spectrum with nonstationarity in latitude (Il4|h From ([3]), the mean can be esti¬ 
mated as (i — T. For each location, we fit a cubic polynomial smoothing spline T min 


from A £JLi{T (L m , i ni t k ) - T m , n (4)} 2 + (1 -A )T.^ 


A = 0.01 since the climate is slowly varying (jCastruccio and Gentonl . 


d/c^ 


-(4) 


with mild penalty term 


20161 ). and we denote 


by T = (T ljl; .... T M>N ). To generate a simulation, the following steps are required: 


• generate e Lm (c) 4 V(0,1 - ¥T m (c) 2 ) with <p Lm (c) as in ([HD, 

• compute H im (c) with (fT3l) and (1T4T) . 


• compute H r (L m ,£ n ,tk) with (ED, 

• compute £ r with (J4J) , 

• obtain the surrogate run as T + e r . 


Once the parameters have been estimated, a common laptop can generate hundreds of 
surrogate runs almost instantaneously with the aforementioned steps. 

In Figure [5] we show a comparison of the six runs from the climate model ensemble and 
the surrogate runs in terms of temperature series near London. In panel (a), the six climate 
model runs are compared with six surrogate runs (offset by 2 C° to avoid superimposition). 
The two groups show the same trend and the same variance, but the statistical model allows 
to generate more runs, so that it is possible to have a better assessment of the temperature 
uncertainty at a given year. (In this context, by “uncertainty” we mean “uncertainty due to 
initial conditions”. We do not consider the uncertainty due to physical parameter calibration 
or forcing scenario.) In panel (b), we see how having just six climate model runs is poorly 
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Figure 5: Comparison of climate model output with surrogate runs, (a) The six realizations 
of the climate model near London (in gray) are shown against six surrogate runs (in red, 
offset by 2 C°). (b) Histogram of the distribution of temperature for the year 2020 for the 
100 surrogate runs. The gray crosses above represent the six realizations from the ensemble 
for the same year. 


informative of the projection uncertainty for 2020, while with 100 surrogate runs it is possible 
to have a better assessment. 

Although a comparison on a single location does not inform about the ability of the 
statistical model to capture the spatial variability, it is possible to produce animations of 


surrogate runs to detect if the spatial patterns are qualitatively consistent. 


Genton et al. 


(120151) have discussed in detail how climate model output and statistical surrogates can be 


compared in the case of three dimensional annual temperatures by using a virtual reality 
environment. In this work, we produce movies for one climate model run and a surrogate run 
(both in the supplementary material), which qualitatively shows similar large-scale features. 
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6 Conclusion and discussion 


In this work we introduced a new class of spectral models that is able to incorporate geo¬ 
graphical information to capture the nonstationary behavior of global data across longitude. 
We further introduced a nonstationary structure across latitude that allows for a more flexible 
and general description of the dependence among different bands. The evolutionary spectrum 
model we developed vastly outperforms axially symmetric models, showing improved perfor¬ 
mance under common model selection metrics and the estimation of the contrast variances. 
By using appropriate diagnostics, we show how this model is able to capture patterns across 
longitude that would be constant under the assumption of axial symmetry. The proposed 
model can be also used to incorporate further geographical information, such as orography, 
or can be applied to other physical quantities whose dynamics are known to be influenced 
by large scale geographical features, such as precipitation or winds. 

The likelihood of the proposed model can be written in a computationally convenient 
form, which is almost as fast as in the axially symmetric case and can be successively ap¬ 
proximated with a highly parallclizable algorithm while still preserving the main space-time 
structure, as shown in the diagnostics. While in this work the approximation blocks and 
the order of approximation (time, longitude, latitude) have been suggested by the particular 
problem, the multi-step approximation presented can be applied to any large space-time 
data set where the nature of the problem suggests blocks: for example, in functional Mag¬ 
netic Resonance Imaging the brain can be naturally divided into regions of interests when 
monitoring cognitive tasks. The analysis was performed on a state-of-the-art workstation, 
allowing distributed computing to optimize the efficiency, and achieving a fit. in less than one 
day for more than 20 million data points. This model consists of 1234 spatial parameters 
and 121824 temporal parameters (three for each location), thus achieving a compression rate 
of 3:100, which is vastly superior to traditional compression algorithms which can achieve 
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at most a 1:5 rate. The fit requires substantial computational power, but the estimates can 
then be used to generate surrogate runs almost instantaneously on a laptop. 

Estimating all the parameters at once would require maximizing a likelihood over more 
than 100,000 parameters for more than 20 million data points, a extremely challenging task 
to perform within a reasonable time even with the most advanced computational facilities. 
Thus, we devised a step-wise estimation procedure with plug-in estimates from previous 
stages, which result in error propagation across stages that needs to be detected and miti¬ 
gated. Bias propagation can be detected with diagnostic figures such as [3] and SI and can 
be mitigated with an intermediate 2-band step before the latitudinal modeling to adjust the 
single band point estimates. Estimation uncertainty propagation in this context is of less 
concern, as given the considerable size of the data set, the estimated standard deviation is 
several orders of magnitude smaller than point estimates and the bias largely dominates the 
error propagation. 

Despite the substantial improvements in flexibility, statistics-based compression is intrin¬ 
sically dependent on the statistical model assumptions. The proposed methodology cannot 
generate surrogate runs that substitute the climate model, as the complex nonlinear dynam¬ 
ics of annual surface temperatures cannot be fully represented by a Gaussian process. As for 
emulators, our statistical model is to be regarded as a useful stochastic approximation that 
could help climate model users to test initial scientific hypotheses, but should not be used 
to perform a full geophysical investigation. 
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